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We study pattern formation on the plane transverse to propagation direction, in a 
ring cavity filled with a Kerr-like medium, subject to an elliptically polarized incom¬ 
ing field, by means of two coupled Lugiato-Lefever equations. We consider a wide 
range of possible values for the coupling parameter between different polarizations, 

B, as may happen in composite materials. Positive and also negative refraction index 
materials are considered. Examples of marginal instability diagrams are shown. It is 
shown that, within the model, instabilities cannot be of codimension higher than 3. 

A method for finding parameters for which codimension 2 or 3 takes place is given. 

The method allows us to choose parameters for which unstable wavenumbers fulhll 
different relations. Numerical integration results where different instabilities coexist 
and compete are shown. 
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PACS: 05.45.-a, 42.65.Hw, 42.70.Mp 


INTRODUCTION 


Spatiotemporal patterns in non linear optical systems, along the plane transverse to light 
propagation, have been widely stndied both theoretically and experimentally [ll, 1^. Stndies 
of optical patterns have some common features with the analysis of pattern formation in 
other physical systems, but there are also some specihc aspects, such as the role of diffraction 
and the vectorial degree of freedom associated with light polarization. 

Patterns taking into account the vectorial degree of freedom of incident helds were an- 

n 

alyzed in [3[ for isotropic positive refractive index material (PRM) with third order non- 


linearities, i.e., a Kerr medium, and mainly for a specihc value of the nonlinear 


parameter 


{B = 1.5). In composite materials an enhancement of nonlinear polarizability jd], and also 
a wider variety of nonlinear parameters may be expected. 

An example of composite materials are negative refraction index materials (NRM); they 
are materials with periodic inclusions which allow the experimental observation of novel 
optical properties, and for which several applications have been proposed Q|. For standard 
PRM, there are arguments to neglect magnetic response, but these arguments do not hold 
for NRM [^. It has been also shown that an NRM can develop a macroscopic effective 
nonlinear magnetic response js]. Negative diffraction is also expected in NRM, but this 
property can also be obtained in regular, periodic refractive index materials j^, Q . Soliton 
formation under zero or negative diffraction has already been studied [h|. 

Here, we extend the study of pattern formation in a ring cavity under arbitrary polarized 
beds, so that it is valid for composite materials, either PRM or NRM, where other values 
of the nonlinear parameter (different form B = 1.5) may be expected, and where nonlinear 
magnetic response may or may not be present. 

We present a method for the analysis of eigenvalues that allows the derivation of some 
exact and general results. The method allows us to hnd parameters for codimension 2 Turing- 
Turing (where two different transverse wavenumbers destabilize simultaneously), Turing- 
Hopf and codimension 3 Turing-Hopf-Turing instability. Also, following our analysis, it can 
be shown that codimension higher than 3 is not possible within the model. 


These results are then used in order to numerically integrate equations, and some results 
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FIG. 1: Scheme of possible ring cavities. In the hrst one, the nonlinear material has length 
L. In the other one, the material hlls the cavity and L is the roundtrip length. Both are 

described by the same equations. 

that are not found in codimension 1 situations are found. 


II. THE SYSTEM 


The system under study is essentially the same as in 12|] with the addition of the trans¬ 
verse spatial dependence. We consider a ring cavity with plane mirrors hlled by an isotropic 
material or metamaterial with a third-order Kerr-like nonlinear response. Two possible 
sketches of the system are shown in Fig. [1] 

The held inside the cavity is described by a plane wave of arbitrary polarization, modu¬ 
lated by a slowly varying envelope. We assume that the electric and magnetic helds are in 
the x-y plane and the wave propagates in the z axis. We study the cavity close to resonance. 

Based on the work by Zharov ef a/ js], we allow the material to have a nonlinear magne¬ 
tization, which depends on the magnetic held. 

Light propagation in a Kerr-type PRM can be described by a nonlinear Schrodinger 
equation, and the same equation can be extended to NRM [l3|. This equation can be used 
to obtain the behavior inside a cavity. Taking into account the magnetic response, and 
applying the same process, we obtain four nonlinear Schrodinger equations (two for the 
envelopes of electric helds and two for the envelopes of magnetic helds, dehned in the plane 
perpendicular to the z axis). It can be shown that the magnetic held remains proportional 
to the electric held. So, the system is well described knowing only the electric held. The 
procedure is analogous to that performed in 
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III. EQUATIONS 


After a change of variables, two coupled Lugiato-Lefever [15| equations, describing the 
left and right circularly polarized held amplitudes inside the cavity, can be obtained: 

dA± 


la 


dt 

i^±r 


— Ain± — (1 4“ + ihV_|_Aj- + 

B\ , . / B' 


1 + 


A, 


( 1 ) 


2 J V 2 

where all cuantities are adimensional, time and transverse coordinates have been normalized; 
A± are the normalized amplitudes of the electric held with circular polarization (see 12|]), 0 


is related to the cavity detuning; a as the sign of XmV‘^~^X^e ■’ with Xe]m being the transforms 
of xxxx component of the third order nonlinear electric and magnetic tensors evaluated at 
~<^o) and rj the inverse of the impedance. The transverse Laplacian, refers to 
the second derivatives with respect to the adimensional coordinates x' = x/l and y' = y/l, 
where / is a characteristic distance (see [l^), h = ±1 is the sign of the dihraction ehects, 
which, in our model, is the same as the sign of the refractive index. Notice, however, that 
negative refractive materials are not necessary for negative dihraction: negative (and zero) 
dihraction resonators can be obtained in negative (or zero) ehective length cavities built by 
means of curved mirrors, see j^. Chapter 6], or by means of a spatially modulated refractive 
index material, see 101. 

The nonlinear parameter B is related to components of the polarization and magneti¬ 
zation tensors that measure the coupling between orthogonal polarization (the nonlinear 
parameter for the electric case is dehned in [^, and the generalization for magnetic non- 
linearities is explained in 12|). Theoretical models predict (see la, p. 227]): = 3/2 in 

materials where nonlinearity is due to molecular orientation ehects; and 5 = 2/3 for elec¬ 
tronic response far from resonance. However, in experiments with Si02 subject to relatively 


long pulses, a value of B as low as 0.244 was measured 


17l |. which was explained as the 


ehect of the competition between electronic and nuclear nonlinearities. Also, the inclusion 
of small spherical particles inside a material, one or both having third order nonlinear re¬ 


sponse, would result in a material w. 
B may take a large range of values 


lere nonlinear ehects might be greatly enhanced, and 
The inclusion of magnetic nonlinear ehects in the 
analysis gives more hexibility to the possible values for B. 

In general, we have a = 1. The less frequent case of a = —1 is equivalent to a = 1, and 
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0 and h with reversed signs, as can be seen by taking the complex conjngate of Eq. o. 
When we have only electric nonlinearities, the case a = — 1 corresponds to a self-defocusing 
material. In the following, we assume a = 1. We also assume that |0| < \/3; within 
this choice, bistable symmetric solutions are not present and changes in 0 do not modify 
qualitatively the results. 

In the rest of our work, we will limit our numerical results to the case 0 < i? <2, and 


where x^m same sign as x^e ■ 

Eq. ([1]) is robust in the sense that not exactly matching impedances can be allowec 


and small dissipation can be taken into account if the normalization is changed, see 


12 


Also, it can be seen that the equation may be still valid for greater values of the detuning 
(new terms can be treated as losses), and diffraction in the linear medium can be taken into 
account redehning the transverse coordinates x and y. 


IV. HOMOGENEOUS SOLUTIONS AND STABILITY ANALYSIS 


Possible homogeneous solutions of Eq. ([T]) were analyzed in 12|, where a classihcation in 
terms of the number of saddle-node and pitchfork bifurcations was presented. They can be 
found solving 


I in 


(1 - = 



( 2 ) 


Where we have dehned the homogeneous solution intensities of the left and right circularly 
polarized components as I± = (where As± are the stationary homogeneous solutions of 

Eq. ([I])), the input intensity as and the polarization as 0 = \Ain+\^/hn 

(the polarization 0 is related to the ellipticity x by 0 = cos^(x/2)). 

For linearly polarized input held (0 = 1/2), there is always a symmetric linearly polarized 
solution, for which = J_. Also, a pitchfork bifurcation may take place at = I' 
producing an elliptically polarized asymmetric solution, where the upper and lower branches 
correspond to either or J_. This new solution may end at = /" > I' (this happens 
if 50 > 2-\/l — B and 5 < 1 and is exemplihed in hg. [2l upper row, for B = 0.9) or may 
not end, i.e. /" = cxo (this happens if 5 > 1 and is exemplihed in hg. [2l upper row, for 
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B = 1.5). 

Instead of symmetric and asymmetric solntions, for elliptic polarization we have contin- 
nous and discontinnous solutions. Continuous solution is present for any value of I in, and 
a discontinuous solution may appear at a given value of the input intensity. A polarization 
0 > 1/2 favors the right circular component for the continuous solution; but the discon¬ 
tinuous solution behaves against intuition, since for 0 > 1/2 we have that /_ > /+. The 
discontinuous solution suddenly starts at like in £g. |3l upper row. It may disappear 

at a second value lin = I". Depending on the parameters, there are three possible situations: 
discontinuous solution absent (for example, for I' —)■ oo), bounded (/' < both hnite, like 
in[3l upper row for B = 1.2), or left unbounded (/' hnite and I” oo, £g. [3l upper row, 
B = 1.5). 

Some basic features of the homogeneous solutions can be analyzed by considering the 
evolution of the perturbations 'ijj± dehned as 


A± = As± 0± 


( 3 ) 


Replacing in ([T]), linearizing and taking the Fourier transform (on transverse coordinates), 
we get 

f Re('0_|_ -f '0_) \ f Re(-0+ -|- 0-) ^ 

lm(0+ '0_) 


dt 


lm('0+ -f 0_) 
Re(0+ - 0_) 

Vlm(0+ -0-) / 

with the linear matrix L given by 

/ —1 9k — aS 


= L 


Re(0+ - 0_) 

V lm(0+ - 0_) / 


( 4 ) 


L = 


( 5 ) 


0 aDB/2 \ 

3aS-9k -1 aD{2-B/2) 0 

0 aDB/2 —1 6k —aS 
-3aD 0 a{l-B)-ek -1 

where S' = J_|_ -I- D = 6*^ = 0 -|- hk'^ and k is the wavenumber of the perturbation 


V 


/ 


the sign of the non linear term (called r; in 
detuning. 

The eigenvalues have the form 


Matrix L has a similar form to the one derived in j^, Eq. (13); one difference is that here 


does not have to be equal to the sign of the 


A±± = -1 ± \JFi± 


( 6 ) 
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where Fi and F 2 are real second order polynomials in 6k'- 

Fi = —6l + bidk + Cl 

F 2 = a2dl + b 26 k + C2 ( 7 ) 

with 

bi = S{3-B/2) 

Cl = S^{B/2 - 2) + D^B/2{1 - B) 

02 = ^^(1 + B/2f - 4 BD‘^{1 - B/2) ( 8 ) 

62 = -2S^{1 + B/2f + 3BSD^{3 -2B + B‘^/A) 

C2 = {S^ + + S^D^B{AB - 5 - 5 ^) 

The homogeneous steady state solution becomes unstable when the real part of one 
of the eigenvalues becomes positive. These instabilities are analyzed numerically in the 
next section, for linear and elliptically polarized input helds. In Sect. IVII we present an 
analytical approach to exactly determine the values of the parameters for specihc situations 
(codimension 2 and 3). 


V. OVERVIEW OF INSTABILITY REGIONS 

In this section we present a general picture of possible patterns and instabilities that can 
occur for different values of the parameters. The parameters are a, the detuning 0, the 
sign of the refraction index h = ±1, the non linear parameter B, and the polarization 0. 
We take a = 1 and 0 = 1 (different values of 0, in the range |0| < \/3 do not produce 
qualitatively different results). Both values of h can be represented in the same marginal 
stability diagram (note that the value of h is not relevant for the shape of the homogeneous 
solutions). 


A. Linear polarization 

From the stability analysis of the symmetric homogeneous solution we obtain the marginal 
stability curves shown in Fig. [2] middle row. In each case two unstable tongues appear, both 
of them are Turing type instabilities. The lowest value of included in the left tongue 


B = 0.9 


B = 1.5 



hn hn 





FIG. 2: Upper row: Homogeneous solutions, /+ (and /_) against Un- Thin curves 
correspond to regions where solutions become unstable under homogeneous perturbations. 
Thick curves correspond to solutions that are stable under homogeneous perturbations. 
Middle and lower row: marginal stability curves, against hk‘^ for symmetric (middle 
row) and asymmetric (lower row) solutions. Asymmetric solutions do not exist in the 
striped region. In all cases, 0 = 1/2. Left column: B = 0.9; Right column: B = 1.5. 


diminishes as B increases, while the right tongue does not depend on B (see Fig. O middle 
row). The point where the left tongue crosses the value fc = 0 corresponds to lin = // 
i.e., it is the point where the symmetric solution becomes unstable under homogeneous 
perturbations, and the pitchfork bifurcation takes place. It is known that, for h = 1, for 
values of close and above the instability threshold of the right tongue, an hexagonal 


pattern appears 


18| . Further increase of the input intensity gives place to oscillating 


hexagons, quasiperiodicity and optical turbulence jl9 |. 


For an NRM {h = —1), close to the instability threshold of the left tongue, a labyrinthic 
pattern is formed at short times when starting from random initial conditions (see |3|). For 
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large times, the system evolves to the homogeneous asymmetric solution. A competition 
between two regions takes place, one with /+ > /_ and the other with J_ > /_|_ (this case is 
illustrated in Fig. 4 of Ref. ^). The marginal stability curves for the asymmetric solution 
are shown in Fig. |2] lower row. The asymmetric solution is always unstable for h = 1. For 
h = —1, there is a range of values of An for which it can be stable. 


B. Elliptic polarization 


B = 1.2 B = 1.5 






FIG. 3: Homogeneous solutions, /+ and /_ against An (upper row) and marginal stability 
curves. An against hk"^ for contiuous (middle row) and discontinuous (lower row) solutions. 
Discontinuous solutions do not exist in the striped region. Labels “H” and “T” refer to 
Hopf or Turing instabilities. In all cases, 0 = 0.6. Left column: B = 1.2; Right column: 

B = 1.5. 

As the ellipticity is increased, starting from 0 = 1/2, for h = —1 the left tongue of 
Fig. [2] middle row is transformed into a closed bounded region whose size decreases until 
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disappearing. For 0 = 0.6, the continuous solution is always stable for h = —1 and for any 
value of B, as can be seen in the marginal stability curves of Fig. [3] middle row. It can be 
shown that, as B is increased, the right tongue changes its shape and is transformed into 
three tongues. Two of them correspond to Turing type instabilities and the central tongue 
is related to oscillatory in time and usually periodic in space Hopf instability (also known 
as Jo type in the notation of Cross and Hohenberg [^h 

The stability analysis of the discontinuous solution shows that, again, for h = —1 and 
0 = 0.6, it is always stable. For h = 1, this solution is always unstable for some k. See Fig. 
|3] lower row. 

In general, similar plots are obtained for other values of (j). Nevertheless, a more detailed 
analysis of the eigenvalues in the plane determined by F 1 -F 2 in eq. (Ej) allows the derivation 
of more general results and the identihcation of some special cases, as explained in the next 
sections. 


VI. DETERMINATION OF INSTABILITY TONGUES 

Since Fi and F 2 in ([ 6 ]) are real quantities, it can be shown that, if one eigenvalue becomes 
positive, then A++ should be positive. Therefore, in order to study stability, it is enough to 
analyze the sign of A++. The analysis is simpler if, instead of describing the unstable zones 
in 9k-Iin or hk^-Iin diagrams, we hrst look at unstable zones in the F 1 -F 2 plane. 

Using that A++ = —1 + Fi + we can see that, for F 2 < 0, points (Fi, F 2 ) that are 
to the right of the line F 2 = 4Fi — 4 have Re(A++) > 0 and Im(A++) 7 ^ 0, therefore, it is an 
oscillatory unstable region. For F 2 > 0, points to the right of the parabola F 2 = {1 — Fi)^ 
with Fi < 1 (that is, the left branch of the parabola) have Re(A++) > 0 and Im(A++) = 0, 
so this region is stationary unstable. The rest of the plane is stable, see Fig. 01 

We are interested in the possible values (Fi,F 2 ) as 9k changes and other parameters 
are hxed. In a marginal instability diagram, changing 9k represents moving through an 
horizontal line. So, if there is a value of 9k for which (Fi,F 2 ) falls in an unstable region of 
Fig. m then, for that value of 9k, in the marginal instability diagram we will be inside an 
unstable tongue. Since Fi and F 2 are second order polynomials in 6 '^ (Ej), the relation can be 
inverted and F 2 can be written as two functions of Fi: F 2 u,i{Fi), which are properly dehned 
in Appendix A. 
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FIG. 4: Unstable regions in the F 1 -F 2 plane. 



f, hk‘ 

FIG. 5: Left: F 1 -F 2 plane with the unstable regions and the curves F 2 u,i{Fi)] arrows 
indicate the direction of growing 6k, in a. situation where the maximum number of 
intersections is obtained. Right: corresponding marginal stability diagram I^n against hk'^, 
points indicated by a, b, c, d, e and f correspond to the intersection points of the left plot. 

Parameters are: 0 = 1, i? = 1, /*„ = 7 and (j) = 0.8. 

An example of the relationship between marginal instability diagrams and Fi — F 2 plots 
is graphically presented in Fig. [5l In the hgure, the values of the parameters were chosen 
in order to get the maximum number of intersections. The left window shows values of 
Fi and F 2 calculated with ([7]) as dk changes and all other values remain hxed. The right 
window shows the related points in the marginal instability diagram. The hgure shows three 
unstable ranges of 9k, the middle one is oscillatory unstable and the other two are stationary 
unstable. This means that the middle tongue in the marginal stability diagram corresponds 
to a Hopf instability and the others to Turing instabilities. This is a general behavior: it 
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can be shown that there can not be more than two tongues related to a Turing instability 
and one Hopf instability tongue. Also, if there are three tongues, the middle one is the one 
related to Hopf instabilities. 

The intersection points of, for example, the hrst tongue and a horizontal line (which 
represents a constant value of /*„), identihed by ‘a’ and ‘b’ in the right window of Fig. |5l 
get closer as decreases, until they merge in one point at the instability threshold. When 
we are at an instability threshold, the curve F 2 u,i{Fi) is tangent to the border of an unstable 
region. The derivation of all these results is mathematically involved and is sketched in 
Appendix A, where other results (most of them intermediate results) are also derived. 


VII. CODIMENSION 2 AND 3 

From the previous analysis we know that we can have, at most, codimension 3 (Turinng- 
Hopf-Turing), i.e., three modes with different wavenumbers that become unstable for the 
same value of lin- We can also have Turing-Turing codimension 2 and Turing-Hopf codi¬ 
mension 2. Figure [6] shows examples of all possible cases of codimension 2 and 3 in a Fi — F 2 
plot and in its corresponding marginal stability diagram. 

Having in mind quite general mathematical properties and constraints given by the phys¬ 
ical system, we are able to derive parameters for codimension 2 and 3. The values of S', 
D and B (that determine the coefficients of Fi and F 2 in ([7])), for which a codimension 2 
or 3 occurs, do not depend on 0. Since 6k = Q + hk^, a change in 0 produces a shift in 
the marginal stability diagram. We can, in principle, take a value of 0 for which instability 
thresholds under consideration are to the left {h = —1) or to the right {h = 1) of k = 0. 

Conditions that parameters should meet in order to have codimension 2 and 3 are derived 
in Appendix B. We summarize the main results here. 

We call Dtt and Stt the values of S and D for a Turing-Turing codimension 2: a similar 
notation is used for the other cases. For 0.848 < B < — 1.028, we have the case of 

Turing-Turing codimension 2. The values of S and D can be found analitically: 


Stt — 


2^/B{3B/2-l) 

B{l-B/2) 


(9) 
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FIG. 6: Curve F 2 against Fi (left) and marginal stability plot (right). From top to bottom: 
Turing-Turing codimension 2, Turing-Hopf-Turing codimension 3, and Turing-Hopf 

codimension 2. 


Dj'x — 


{1 + B/2)S 

V-BV2 + 3S' 


( 10 ) 


For B = Be, we have codimension 3. Stht and Dtht are given by Stht = Stt{Bc); 
Dtht = Dtt{Bc). 

Parameters for Turing Hopf codimenision 2 situations are harder to determine, see the 
second part of Appendix B. After some algebra, we hnd that conditions for codimension 
2 Turing-Hopf situations are met only if the roots of a given polinomial P{r), which, once 
S and B are hxed, is fourth degree in an auxiliary variable r (related to the difference 
between Fi and its maximum value), has a double real and two complex conjugate roots (or 
two double roots; which only happens for B = B^ and corresponds to the codimension 3 
situation previously described). For every value of S and for B > Be, we can numerically 
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5 

FIG. 7: Number of real roots from of P{r). Regions where there are no possible 
codimension 2 situations (if there is no value for D, or if it is greater than S) are also 
shown. Curves in the Fi — F 2 diargram become tangent and thus, Turing-Hopf 
codimension 2 takes place for parameters given by the thick line of this hgure. Notice that 
no codimension 2 Turing-Hopf can occur for B < Be 

compute the roots of that polynomial. FigureHshows the number of roots in the plane S-B. 
The thick line represents the set of points for which codimension 2 Turing-Hopf situations 
take place. It is worth mentioning that for every allowed B there are two possible values of S\ 
the one to the left (right) of the vertex has a value oiQ + hk^^ greater (lower) than Q + hk‘^ 
(where kn is the expected wavenumber for one of the Hopf instabilities, and similarly for 
kx), so, for h = 1 the wavevector related to Turing instability is smaller (greater) than the 
one related to Hopf instability. Finally, taking a point of this curve the values of S, B and 
D can be determined, see appendix B. 

In all cases, following the derivations shown in the appendices, and choosing a value 
of 0 we get the unstable wavenumbers, for instance kxi and /cr 2 for the Turing-Turing 
codimension 2 situation. Conversely, we can choose the unstable wavenumbers (for instance 
we can make them fulhll a given relationship) and use that information to properly choose 

0 . 

Once we hnd B, S and D, having in mind that I± = and choosing a value for 0, we 
can hnd the input intensity lin and polarization (p for which codimension 2 or 3 takes place 
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in a straightforward way (just replacing all known values in Eq. (|2]) an solving two coupled 
linear equations). 


VIII. NUMERICAL INTEGRATION RESULTS 


Numerical integrations of Lugiato-Lefever equations have been extensively performed in 
previous reports. The novelty here is that we will exploit the results from previous sections 
in order to hnd parameters for codimension 2 and 3 in a straightforward way. The purpose 
of this section is to have a quick look at possible situations that may occur when patterns 
tend to emerge in codimension 2 or 3. 

In codimension 2 Turing-Hopf situations were analyzed for the special case B = 1.5. 
They found out that an hexagon related to a Turing instability dominated at long times, 
although the Hopf instability dominated at short times. Also, in 2l|] they analyzed a Turing 
- Turing codimension 2 instability and found that different patterns related to competition 
of unstable wavenumbers might take place. 


A similar research, for a different system, was performed in 22| , where pattern formation 
situations are analyzed in a Belusov-Zhabotinsky reaction, and codimension 2 Turing-Hopf 
may occur. They found out that in codimension 2 situations, patterns related to both 
instabilities coexist for quite long times, but eventually one dominates. The exception occurs 
in one dimension when destabilizing modes are resonant (for instance, the wavelength of one 
instability is an integer times the wavelength of the other instability), in that case, both 


23| for a reaction-diffusion 


unstable modes may coexist. Similar results where found in 
model where also chaotic situations are allowed. 

In optics, two coupled Kerr-like systems (specihcally, two liquid crystal light valves) 
where analyzed both theoretically and experimentally 24j-l26| . Turing-Hopf codimension 
2 situations were reported. For some parameters, unstable wavevectors where resonant, 
and a far held composed of two octagons (whose radius where the wavevector modules 
of the unstable modes), one of them rotated tt/S degrees respect to the other, was found 
251]. Turing-Turing codimension 2 or higher codimension was not allowed since a linearly 
polarized system was studied (and, from the dynamical point of view, the system was two 
dimensional, i.e. instead of the matrix in ([5]), they had a two by two matrix). 

Here, we are interested in situations where the sum of unstable modes related to one 
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instability may contribute to a mode related to another instability. For instance, if one 
instability is related to hexagonal patterns with some orientation and the second one has 
a wavelength y/S times greater, we expect the second one to form an hexagon y/S times 
larger, and rotated tt/G degrees respect to the first one (so that the sum of wavevectors of 
the hrst instability should contribute to the other instability). Notice that parameters for 
which wavenumbers of the different instabilities fulhll desired relations can be found taking 
into account the calculations performed in previous sections. 

Results of our numerical integrations j^, are the following: At short times, all unstable 
wavevector coexist (so that |A±p is composed of two or three rings in the far held), then 
rings of unstable wavector become thinner and the intensity of one of them becomes much 
greater than the others. After that, different situations may occur. 

In some nonresonant codimension 2 Turing-Turing situations, we found that at long times 
a ring of unstable modes (of radio ku in the far held) with wavevectors diherent but among 
the values of the unstable modes dominated {kxi < kjj < k'r 2 ), and the near held was 
composed of domains of ordered hexagonal patterns (for example, A+ with up hexagons, 
and A_ with down hexagons). For some parameters, instead of a ring of modes, in the far 
held an hexagon (with kxi < ku < kT 2 ) was formed, and an hexagonal pattern arose, with 
a unique orientation, in the whole near held. 

Taking an adequate value of 0, it is possible to make unstable wavenumbers fulhll the 
desired ratio. Simulations with kri = fcr2/\/3 were performed. For some parameters, kT 2 
dominated and a regular dodecagon took place in the far held. Finally, putting an initial 
condition that was the steady solution plus an hexagonal pattern related to the smallest 
wavenumber, it could be seen that both unstable wavenumbers grew, forming an organized 
hexagonal structure in the near and in the far held. The same result was found even when 
the input intensity was slightly lower than critical intensity. See Fig. |8l where, for |A_|_p, the 
coexistence of two unstable wavenumbers can be found even in the near held. The stability 
analysis for this case is shown in Fig|9l it can be seen that there are 2 unstable wavenumbers, 
i.e. two values of k for which i?e(A++) is not negative. 

Codimension 3 situations were also analyzed. Parameters were chosen so that the ratio 
between the greatest Turing wavelength and the smallest one was y/3 {kxi < ku < kj ’2 = 
yf^kTi)- At the initial stage, we found that all unstable wavenumbers where activated (see 
the left plot on hg. [T0|) . forming three concentric circumferences in the far held. The smallest 
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FIG. 8: Left column: |A+P, right column: |A-P, up: near field, down: far field, at t=220, 
for the case of a resonant Turing Turing codimension 2 where initial conditions have an 
hexagonal modulation. The same results are found at larger times. Parameters: h = 1, 
B = 0.94, 0 = 1.499, 0 = 0.886 and = 2.799, t = 220. 



FIG. 9: Real (solid line) and imaginary (dotted line) parts of A++ as a function of k, for 
the parameters of the previous hgure. Vertical lines show unstable wavenumbers (fcn and 
kT 2 ), related to steady perturbations. The ratio among them is a/S. They closely match 
the numerical integration results of previous hgure. 


Turing wavenumber grew faster, and turned into an hexagon. After that, another hexagon. 
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related to the greatest Turing wavenumber also appeared, see Fig. [TO] middle plot. At odds 
with £g. [HI regions where instabilities with different wavenumbers dominate are spatially 
separate in the near field. Finally, a crown of modes got destabilized (right plot on fig. [TOli . 
The stability analysis for this case is shown in Fig[TTl Notice that the intermediate unstable 
wavenumber is related to Hopf instability, i.e. it has /m(A++) 7 ^ 0. 

In another example of resonant Turing-Hopf-Turing Codimension 3, in which kT 2 = 2 /cti 5 
a similar final situation was observed. For nonresonant cases, there were found steady 
situations similar to fig. [lOl middle graph, but where one Turing wavevector (making an 
hexagon or a ring in the far field) prevailed for A+ but the other prevailed for A_. 

For Turing Hopf codimension 2, situations similar to Turing Turing codimension 2 were 
found, both in resonant and nonresonant situations. Also, square patterns took place in 
some numerical integrations. 


IX. CONCLUSIONS 

Taking into account the polarization degree of freedom of light, and having in mind that 
B, which measures the nonlinear coupling between different polarizations, could take a broad 
range of values, we presented a study of instabilities and patterns that might show up in a 
cavity filled with a Kerr-like nonlinear material with positive or negative refractive index. 

A method for finding codimension 2 and 3 situations (where wavenumbers of different 
modulus might destabilize) was shown. It can be found that for 0.848 < B < Turing- 
Turing codimension 2 may occur; for B = B^ Turing-Hopf-Turing codimension 3 may take 
place, and for B > B^. there may be situations for Turing-Hopf codimension 2. Fixing 
only the value of B, the method allows us to find all other values of the parameters for 
codimension 2 or 3. It allows also to see that, for a given intensity, there cannot be more 
that three instability regions in a marginal instability plot (one of which has to be related 
to a Hopf instability), and that codimension higher than 3 cannot occur. 

Since the method allows us to know some parameters with any degree of precision (instead 
of performing a numerical search and changing the parameters until such situation shows 
up), and choose others at will, it is a useful tool in the study of codimension 2 or 3 on the 
model. Specifically, resonant situations, where the ratios between unstable wavevectors are 
chosen, can be found. Also, it might be useful for the underestanding of pattern formation 
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in other systems as long as the linear stability analysis presents eigenvalues with the form 
of eqs. ([6]) and ([7|). 

Numerical integration results show some new situations of pattern coexistence and com¬ 
petition. 
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APPENDIX A: PROPERTIES OF INSTABILITY TONGUES 

We consider that the possible values (^ 1 ,^ 2 ) are those given by these polynomials with 
the free parameter 6k, and with hxed coefficients, i.e., we consider fixed values of S and D, 
that correspond to a hxed value of the input intensity for a determined homogeneous 
solution. From ([7]) we see that Fi takes a maximum value given by Fim — h\/A -\- Cl- 
By dehnition, S > D, so it can be shown that 02 > 0, and that Fi + \/p 2 —00 for 

\6k\ — )■ 00 . This means that, for large \6k\, Re(A++) = —1. Therefore, the range of unstable 
wavenumbers is bounded. The case 02 = 0 occurs only ii S = D, that is, for pure circular 
polarization, but this case can be related to the pure linear polarization case, as has been 
done, for instance in j^. 

From ([7]), we can obtain two solutions for 9k as a function of Fi. Using these solutions 
in the equation for F 2 , we obtain: 

F 2 u,i{Fi) = (02 & 1/2 -|- 62 bi/2 + 02 + 0,2 Cl) — 02 Fi + \X\r, (A.l) 

where X = 02^1 + ^2 and r = a/Tim — Fi > 0; indices ‘u’ and stand for upper and lower 
curves respectively. 

The difference between the upper and lower curves is F 2 u{Fi) — F 2 i{Fi) = |Y|r. So, we 
have two values of F 2 for each Fi as long as Fi < Fim and X 7 ^ 0. For X = 0, the curves 

F 

are two overlapping rays that start from Fi = Fim- We can also see that = T2|X|/r, 

so F 2 u {F 21 ) has negative (positive) curvature. 

To obtain the instability points, we have to look at the intersection of F 2 u,i{Fi) with the 
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unstable regions of Fig. HI Eq. (lA.lj) can be rewritten 


as: 


{F2-d + = X^{F,m - F,), (A. 2 ) 

where d = 02 6i/2 + 62 &i/2 + C 2 + 02 ci (the same equation holds for F 2 U and F 21 ). To obtain 
the intersection with the stationary unstable region, we replace F 2 by (1 — Fi)^ in flA.2ll . 
We get a 4th order polynomial, so there are at most 4 solutions for Fi. The intersection 
with the oscillatory unstable region is obtained replacing F 2 by 4Fi — 4 in flA.2j) : this gives 
two solutions for Fi. To obtain the maximum number of unstable ranges of values of 9k 
we assume that half of the previously mentioned solutions of Fi correspond to a cross from 
stable to unstable region (as 6k is increased) and the other half to a cross from unstable to 
stable region. So, we have, at most, three unstable ranges of 9k that correspond to three 
tongues in the marginal stability diagrams. 


APPENDIX B: DERIVATION OF PARAMETERS VALID FOR CODIMENSION 

2 AND 3 

For Turing-Turing codimension 2, we need both F 2 U and F 21 to be tangent to the border 
of the stationary unstable region, and this happens only if X = 0 ; in this case F 2 U and 
F 21 are straight lines that overlap. To have codimension 3 (Turing-Hopf-Turing) there is a 
further condition: the end point of the rays, {Fim, F 2 {Fim)), must be on the border of the 
oscillatory unstable region given by F 2 = 4Fi — 4. The case of Turing-Hopf codimension 
2 occurs when F 2 against Fi is tangent to the borders of the stationary and oscillatory 
unstable regions, and X ^ 0. See Fig. [6l 

Turing-Turing codimension 2 and Turing-Hopf-Turing codimension 3 

From the condition X = 02^1 -|- 62 = 0 it is straightforward to obtain an expression 
for ffTOj) . From Eq. flA.2p we have that F 2 = d — 02 ^ 1 , and the intersection 

with the border of the stationary unstable region, given by F 2 = (1 — Fi)^, gives a 2 nd 
order polynomial in Fi. To have the line tangent to the parabola, the discriminant of the 
polynomial should be zero. From this conditions we get the critical value Fic = 1 — 02/2 
from which, using Eq. ([7]), we get the two critical values of 9k- From the zero discriminant 
and Eq. (nnii. we get an expression for ([9]). 
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Then, for a given value of B, there is a unique value of Stt and \Dtt\ where we can find 
Turing-Turing codimension 2. There are some restrictions on the possible values of B. First, 
in order to have Stt real, we have that B > 2/3, but there is a more restrictive condition. 
We need that Fic < Fim in order to have a solution tangent to the unstable border that 
actually touches it. It can be shown that the condition B > 0.848 should be satisfied. 
Second, the Hopf instability should appear for greater values of than the Turing-Turing 
instability. 

We define the distance between F 2 {Fim) and the border of the oscillatory unstable region 
as Z = F 2 {Fim) — {4: Fim — 4) . For Turing-Turing codimension 2 we need Z > 0, and 
for codimension 3, we have that Z = 0 since the point {Fim, F 2 {Fim)) should be on the 
oscillatory unstable border. Since Stt and Dtt are functions of B [see Eqs. ([9]) and (ITOl) ]. 
we can obtain Z as a function of B only. It can be shown that the only zero of Z occurs for 
Be = 1.028. Then, using the value B = Be, we can obtain the parameters Stht = STT{Be) 
and Dtht = Dtt{Bc) for codimension 3. For B < Be, Z > 0, so that Turing-Turing 
codimension 2 is allowed. For B > Be, Z < 0: we still have the two stationary instabilities 
that occur simultaneously for a given value of /j„, but this is not Turing-Turing codimension 
2 since the oscillatory instability appears for a smaller value of (As we will see below, in 
that region there are Turing-Hopf instabilities). 

Expressions for 0 -|- hk'^i and 0 -|- hk ^2 can be found solving Fi = 1 — 02 / 2 , and replacing 
S by Stt and D by Dtt- If there is codimension 3, the value of 0 -|- hk^^j can be found from 
0 -I- hk'jj = hi/2. It can be shown that 2hk‘jj = hk'^i + hk'^ 2 - Once the value of 0 is chosen, 
kn, kTi and kT 2 are fixed. Conversely, once two wavenumbers are chosen, 0 is fixed (and so 
is the third wavenumber, if exists). 


Turing-Hopf codimension 2 

For a Turing-Hopf codimension 2 we require the curve F 2 {Fi) to be tangent to both 
borders of the unstable regions, as shown in the lower row of Fig. | 6 l 

Let us hrst consider the contact point with the border of the oscillatory unstable region, 
i.e., between F 2 i{r) (lA.lIl and the line F 2 = 4Fi — 4 = 4{Fim — r^) — 4. The intersections 
are given by a second order polynomial in r. We require the intersection to be only in 
one point, so that the polynomial discriminant should be zero. From this condition, it is 
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possible to obtain three possible expressions for D as a function of S and B. We will call 
them 

Now, we consider the intersection with the stationary unstable region, i.e., between F 2 u{r) 
and F 2 = (1 — -Fi)^ = (1 — Fim + r^)^. We get a fourth order polynomial in r, which will be 
called P{r). We can have 0, 2 or 4 real roots, and we are interested in the cases of a fourfold 
real root, or double real and two complex conjugate roots, in order to have the function F 2 U 
tangent to the unstable border. It can be shown that the kind of roots that we are looking 
for are possible only for one of the expressions of D mentioned in the previous paragraph, 
say Di{S, B). In £g. [7] we plot the roots of P{r). Its coefficients are calculated for S', B and 
D = D,{S,B). 

The value of Q + hk^ can be found solving = Fim — ^q, where tq is the double 

root of P{r). The value of 0 + hkjj can be found in a similar way. 
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FIG. 10: l^+l^ in the near (up) and far field (down) at different times for the case of 
resonant codimension 3. From left to right, t = 333 (all unstable wavenumbers are 
enabled), t = 482 (an hexagonal structure in the far held) and t = 570. The A- component 
has a similar behavior. Parameters for this numerical integration are: h = 1, B = Be, 

0 = 1.15, 0 = 0.78 and hr, = 4.26. 



FIG. 11: Real (solid line) and imaginary (dotted line) parts of A++ as a function of k, for 
the parameters of the previous hgure. Vertical lines show the values of the unstable 
wavenumbers, two of them (fcn and kT 2 ) are related to steady perturbations, and the 
ratios among them is \/3. The other one (kn) is intermediate between them and is related 
to oscillatory instabilities (/m(A++) 7 ^ 0). kxi, kn and /ct 2 closely match the unstable 

wavenumbers in the previous hgure. 























